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1.  INTRODUCTION 

A  multirate  method  for  integrating  ordinary  differential  equations  is 
one  in  which  different  equations  are  integrated  using  different  time  steps. 
The  motivation  for  such  an  approach  is  obvious — the  desire  to  reduce  the 
integration  time.  Consequently,  the  idea  has  been  mentioned  many  times  and 
probably  used  on  an  ad  hoc  basis  in  a  number  of  applications,  particularly 
in  real-time  applications  where  speed  is  critical.  However,  there  has  been 
relatively  little  study  of  the  problems  involved  in  efficient,  automatic, 
multirate  integration.  In  this  paper  we  will  first  look  at  the  types  of 
equations  that  can  be  handled  efficiently  by  multirate  methods.  Next  we 
will  take  a  brief  look  at  the  theory  of  multirate  methods.  We  will  see 
that  there  are  no  complications  in  the  conventional  theory  of  convergence, 
and  that  the  expected  happens.  Then  we  will  look  at  the  problems  of 
automatic  multirate  methods.  Here  we  encounter  a  number  of  difficulties  in 
the  estimation  of  errors  and  the  control  of  stepsizes.  The  latter  occurs 
because  it  may  not  be  possible  to  estimate  the  error  in  a  particular  step 
until  after  a  number  of  other  steps  have  been  taken,  by  which  time  it  may 
be  very  costly  to  "back-up"  the  integration  to  redo  the  earlier  step. 
Three  different  techniques  which  have  been  studied  will  be  discussed. 

2.  APPLICATION  OF  MULTISTEP  METHODS 

One  area  in  which  multirate  methods  have  been  considered  is  simulation 
{1}{2},  although  there  is  little  documentation  on  the  subject.  Some  work 
has  been  done  on  partitioning  systems  into  subsystems  for  treatment  by 
different  methods,  as  in  {3}  where  the  nonlinear  part  is  integrated  by  a 
conventional  Runge-Kutta  method  and  the  linear  part  is  handled  by  methods 
that  use  precomputed  matrix  exponentials.  In  our  case,  we  are  interested 
in  partitioning  based  on  the  relative  speeds  of  subsystems.  In  such 
applications,  the  engineer  is  often  able  to  use  a  large  amount  of 
"engineering  intuition"  in  the  design  of  an  integrator.  For  example, 
consider  the  simulation  of  an  aircraft.  It  consists  of  a  number  of 
subsystems.   For  simplicity  we  will  consider  two  shown  in  fig.  1. 
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Fig.  1.   Simplified  simulation  model 

The  control  subsystem  reacts  very  rapidly  (being  electronic),  whereas  the 
flight  dynamic  subsystem  reacts  relatively  slowly,  being  a  mechanical 
change.  Consequently,  a  simulation  system  tailored  to  this  model  might 
choose  to  integrate  the  control  with  small  stepsizes,  and  the  dynamics  with 
large  stepsizes.  The  rationale  sometimes  heard  for  this  is  that  "because 
the  dynamics  are  slow,  we  can  break  the  feedback  loop  from  the  dynamics  to 
the  control  and  handle  each  subsystem  separately."  "Breaking  the  feedback 
loop"  refers  to  assuming  that  there  is  no  change  in  the  output  of  the  slow 
system  over  one  of  its  steps  while  the  fast  system  is  being  integrated  over 
several  smaller  steps.  What  is  actually  happening?  Suppose  the  outputs  of 
the  fast  and  slow  systems  are  designated  as  y  and  z  respectively,  and  they 
satisfy  the  differential  equations 


y'  =  f(y,z,t) 

z'  -  g(y»z,t) 


(l) 


(These  could  be  two  systems  of  equations  of  different  dimensions.)  Let  us 
suppose  that  the  stepsizes  used  for  the  y  and  z  integrations  are  h  and  Mh 
respectively,  where  M  is  an  integer,  and  that  the  values  of  y  and  z  are 
known  for  time  value  tw^,  where  t^  =  t*  4-  kh.  The  process  consists  of 
integrating  y  over  M  small  steps  of  size  h  while  z  is  held  constant  at  its 
value  of  z (£«_)•  Then  z  is  integrated  over  one  step  of  size  Mh  so  that  we 
have  information  at  time  t„, n+i\*  Fig.  2  shows  a  possible  response  of  the 
numerical  system  where  f(y,z,t)  =  -100 (y-z).  The  y  values  change  rapidly 
after  z  has  been  updated,  then  settle  to  a  new  "stable"  state  based  on  the 
value  of  z  at  the  beginning  of  the  group  of  M  small  steps.  In  fact,  we 
would  expect  that  the  true  value  of  y  at  the  end  of  the  set  of  M  small 
steps  should  be  determined  by  the  value  of  z  at  the  end,  and  this  suggests 
that  it  might  be  better  to  perform  the  integration  of  z  over  its  big  step 
before  doing  the  M  steps  of  y.  This  would  correspond  to  "breaking  the 
feedback  loop"  between  the  control  and  dynamic  subsystems;  intuitively  less 
appealing  but  possibly  more  accurate.   At  first  glance  it  might  appear  that 
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this  simply  transfers  the  problem  to  z  which  now  will  be  integrated  over 
one  step  Mh  based  on  the  value  of  y  at  the  start  of  that  step.  This  is  not 
a  problem  if  the  integration  method  used  for  z  is  explicit  because  the 
value  of  g(y,z,t)  is  needed  only  from  the  beginning  of  the  step,  but  there 
is  a  more  important  reason  why  it  is  not  a  problem.  If  the  behavior  of  y 
is  fast  but  of  z  is  slow,  there  can  be  very  little  coupling  from  y  to  z. 
That  is  to  say,  the  dependence  of  g(y,z,t)  on  y  must  be  very  small.  Hence 
the  link  from  y  to  z  is  a  natural  place  to  "break  the  feedback  loop."  This 
point  will  be  important  in  later  discussions. 


Fig.  2.   Integration  after  "Breaking  the  Feedback  Loop" 


Wherever  the  feedback  loop  is  broken,  a  constant  value  of  one  variable 
is  being  used  in  the  integration  of  the  other  variable  over  an  interval  of 
length  Mh.   This  "approximation"  by  a  constant  function  will  introduce  an 

error  of  size  O(Mh),   and  its  integration  over  an  interval  of  size  Mh 

2 
contributes  a  local  error  of  O(Mh)  which  causes  a  global  error  of  size 

O(Mh).   The  natural  thing  to  do  is  to  interpolate  or  extrapolate  the  value 

of  z  (or  y)  while  integrating  y  (or  z)   thus  reducing  the  approximation 

error.   If  the  extrapolation  is  done  by  a  polynomial  of  degree  p-1  it 

introduces  an  extrapolation  error  of  0(Mh)P.   In  that  case  the  global  error 

introduced  will  also  be  0(Mh)P,  which  would  be  a  natural  ally  of  an  order  p 
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integration  method.  Section  4  will  sketch  the  proof  of  the  result  being 
suggested  here,  namely  that  if  stable  order  p  methods  are  used  for 
integration  and  order  p-1  methods  are  used  for  extrapolation  or 
interpolation,  the  method  is  order  p  convergent. 

Another  area  in  which  multirate  methods  are  being  considered,  although 
under  another  name,  is  in  evolutionary  partial  differential  equations 
(PDEs)  where  it  is  thought  of  as  a  form  of  local  mesh  refinement  in  time. 
Such  methods  will  probably  be  necessary  if  large  two-  and  three-dimensional 
time-dependent  PDEs  are  to  be  integrated  in  the  presence  of  shock-like 
fronts  because  the  size  of  the  time  step  must  be  sharply  restricted  in  the 
front  but  not  in  regions  of  low  activity.  In  this  paper  we  consider  the 
general  problem  of  systems  of  ordinary  differential  equations  which  could 
arise  from  a  discrete  system  being  simulated  or  by  the  spatial 
discretization  of  PDEs.  In  the  latter  case,  we  should  recognize  that  the 
system  of  equations  may  change  from  time  to  time  as  the  spatial  mesh  is 
adjusted  to  adapt  to  the  needs  of  the  problem,  but  these  particular 
questions  will  not  be  considered  here. 

3.   SPEED  CONSIDERATIONS 

Multirate  methods  appear  attractive  because  they  should  use  less 
computer  time.  The  reasoning  is  that  the  total  time  is  roughly 
proportional  to  the  number  of  integration  steps  taken  for  each  equation. 
Since  a  conventional  method  must  use  a  number  of  steps  equal  to  the  number 
of  steps  for  the  fastest  component  times  the  number  of  equations,  a 
multirate  method  can  reduce  the  number  of  steps  for  the  slow  components. 
However,  multirate  programs  can  be  more  complex  than  straightforward 
methods,  so  we  must  be  certain  that  increased  program  overhead  does  not 
override  any  savings.  The  computational  cost  of  a  numerical  integration 
program  is  due  to  several  factors.  They  are: 

Evaluations  of  derivatives 

Application  of  integration  formulas 

Interpolations/extrapolations  if  a  multirate  method  is  used 

Estimating  errors 

Logic  for  automatic  step  control 
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Repeated  steps  when  a  failure  occurs  and  completed  steps  must 

be  discarded 
Solution  of  implicit  equations  if  they  are  used 

In  this  section  we  will  examine  the  relation  between  the  first  three  of 
these  costs  for  the  pair  of  eq.  (1)  when  they  are  integrated  by  fixed- 
stepsize  methods. 

Let  us  assume  that  the  same  method  is  used  for  each  equation  over 
stepsizes  of  h  and  Mh,  and  that  it  requires  q  function  e  fluations  per 
step.  We  will  consider  the  cost  of  one  compound  step  of  size  Mh.  Suppose 
that  C£,  c  ,  c  ,  and  c.  are  the  costs  of  an  evaluation  of  f,  an  evaluation 
of  g,  an  extrapolation  of  z,  and  an  integration  step,  respectively.  Then 
the  cost  of  one  compound  step  is 


cm  -  <qcf  +  cx>M  +  icg  +  <M  +  l>ci  (2) 

This  assumes  that  only  one  interpolation/extrapolation  is  needed  for  a  z 
value  in  each  integration  step  for  y,  as  would  be  the  case  in  a  multistep 
method.  If  intermediate  values  are  used,  as  in  a  Runge  Kutta  method,  the 
cost  must  be  increased  by  M(q-l)c  .  This  must  be  compared  with  the  cost  of 
using  stepsize  h  for  both  components  over  M  steps,  which  is 


C  -  (Cf  +  cg)Mq  +  2Mc± 
because  no  extrapolations  will  be  necessary.   The  difference  is 


C  -CM«  (qcg  +  Ci)(M-  1)  -M[q]cx  ^ 

where  the  [q]  term  appears  if  a  Runge-Kutta-like  method  is  used.  In 
Runge-Kutta  methods,  c^  is  small  and  q  is  large  (4  to  13,  depending  on  the 
order  of  the  method).   In  this  case  we  can  approximate  the  difference  by 


C  "CM*  <cg  "  cx)^  "  »*   ~   <cx  (5) 

Therefore,  there  is  nothing  to  be  gained  from  multirate  Runge-Kutta  methods 
unless  c  >  c  in  which  case  there  is  a  saving  if 

6      X 
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M  >  1  +  -z *—=- 

Cg  "  Cx  (6) 

If  c   is  only  a  little  larger  than  c  ,  the  slight  savings  for  large  M  will 
g  x 

be  lost  to  the  undoubtedly  higher  overhead  of  a  more  complex  program.  On 
the  other  hand,  for  multistep  methods,  q  is  small,  typically  2.  Also,  the 
number  of  prior  points  used  for  extrapolation  is  naturally  the  number  of 
points  used  in  the  multistep  method  for  two  reasons.  The  first  is  that  we 
are  saving  values  at  those  points  anyway.  The  second  is  that  the  predictor 
step  of  a  multistep  method  is  an  extrapolation;  we  have  chosen  the  number 
of  values  and  the  stepsize  in  the  multistep  method  such  that  the  error  of 
the  prediction  is  not  too  large,  and  this  implies  that  the  error  in  the 
extrapolation  will  also  be  appropriately  small  if  the  same  number  of  past 
values  are  used.  Consequently,  the  cost  of  a  multistep  integration  step  is 
only  slightly  higher  than  the  cost  of  an  extrapolation  because  it  involves 
a  prediction  step  followed  by  a  correction.   Therefore, 


cx  <  ci 


Using  this  in  eq.  (A),  we  get 


C  -  CM  >  q(M  -  l)cg  -  cx  (7) 

Therefore,  we  may  see  a  saving  as  long  as 


M  >  1  +  x 


qcg  (8) 

Eqs.  (6)  and  (8)  indicate  that  multistep  methods  are  more  likely  to  yield  a 
savings   in  multirate  schemes  than  Runge-Kutta  methods,  but  that  unless  c 


g 
hat   thei 

must  be  a  very  large  disparity  in  behavior  of  the  two  components. 


is  large  compared  to  c  ,  a  large  M  will  be  needed,  which  means   that   there 


Similar  comments  apply  to  larger  systems,  but  sparsity  can  have  a 
positive  effect.  Suppose  that  y  and  z  above  each  represent  systems  of  r 
equations,  but  that  the  evaluation  of  f  requires  the  values  of  only  s  of 
the  z  variables.  Then  the  extrapolation  cost  is  scx  whereas  the  cost  of 
evaluation  of  all  of  the  components  of  g  is  re  if  cg  is  the  cost  for  each 
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component.  Thus  the  sparsity  ratio  s/r  multiplies  c  in  eqs.  (6)  and  (8), 
reducing  the  value  M  for  which  savings  can  be  achieved.  Clearly  this 
analysis  can  be  extended  to  cases  in  which  more  than  two  different 
stepsizes  are  employed  and  the  same  general  result  will  be  obtained,  namely 
that  the  ratio  of  the  cost  of  extrapolation  to  the  cost  of  function 
evaluation  is  critical.  Also  note  that  it  is  very  important  where  possible 
different  components  are  evaluated  at  the  same  point.  If,  for  example,  the 
two  components  in  eq.  (1)  were  integrated  with  the  same  stepsize  on 
different  meshes,  say  tQ  +  nh  and  t«  +  (n+l/2)h,  two  unnecessary 
extrapolations  would  be  done  in  each  step,  and  the  method  would  take  longer 
than  if  the  equations  were  "kept  in  step." 

4.   FIXED-STEP  CONVERGENCE  THEORY 

It  is  natural  to  wonder  if  the  use  of  multirate  techniques  can  cause 
problems  with  stability  and  convergence.  This  section  will  show  that  they 
do  not.  Indeed,  the  expected  happens:  if  the  methods  used  for  each 
equation  are  individually  stable  and  p-th  order  accurate,  if  the 
functions — f  and  g  in  eq  (1) — are  Lipschitz  continuous  in  all  variables, 
and  if  the  extrapolation  scheme  is  (p-l)-th  order  accurate — in  the  sense 
that  the  error  in  extrapolation  is  O(hP),  then  the  solution  is  p-th  order 
convergent.  This  result  is  true  for  any  number  of  equations,  each  with  a 
separate  integrator  using  a  separate  stepsize.  We  will  sketch  the  proof 
for  two  equations.  It  is  a  straightforward  extension  of  standard  proofs, 
and  this  sketch  will  make  it  clear  how  the  result  can  be  obtained  in  any 
particular  application  without  obscuring  the  principles  with  generality. 

Consider  the  numerical  integration  of  a  system  of  equations  y'  = 
f(y,t)  all  using  same  stepsize  h  for  the  n-th  step.  It  has  been  shown  in 
{4}  that  the  global  error  e^  at  the  n-th  step  of  a  method  satisfies  the 
equation 

n+l        vnnn'nn  /n\ 

where  S^  is  the  stability  matrix  for  one  step  of  the  method  with  stepsize 
h^  applied  to  the  differential  equation  y'  =  0,  S^  is  the  contribution  of 
terms  dependent  on  the  differential  equation  (and  depends  on  the  partials 
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f  ),  and  cL  is  the   "local  truncation  error."   In   this   case   the   local 
y' '      n 

truncation  error  is  defined  quite  simply;  it  is  the  amount  by  which  the 
true  solution  fails  to  satisfy  the  method.  In  the  case  of  a  multirate 
method,  eq.  (9)  is  applicable  to  each  set  of  differential  equations  being 
integrated  with  the  same  stepsize,  but  d?  is  composed  of  three  components: 
the  truncation  error  due  to  the  method  applied  to  this  set  of  differential 
equations,  the  interpolation/extrapolation  errors  when  other  components  are 
approximated  for  the  evaluation  of  f,  and  the  global  errors  in  the  other 
components  used  to  evaluate  f.  For  the  case  of  the  two  equations  (1),  we 
will  name  these  components  d^,  iz,  and  ez,  respectively.  The  latter  two 
terms  enter  into  eq.  (9)  through  the  function  f(y,z,t),  and  are  always 
multiplied  by  h  f  ,  hence  the  revised  form  of  eq.  (9)  is 


ej+i    -    (S£  +  V^ey  +  h  f   e*   +  hf   i*   +  d£ 
n+1  n  nnn  nzn  nzn  n 


(10) 


A  similar  type  of  equation  will  apply  to  the  global  errors  in  z,   although 

it  will  apply   to  a  different  set  of  mesh  points.   However,  we  can  define 

values  of  the  e^  and  ez  on  the  union  of  all  mesh  points  used  in  any 

equation  to  be  the  values  of  the  global  errors  at  the  last  point  to  which 

the  component  was  actually  integrated.   Then  a  form  of  eq.  (10)  applies   to 

all  components  at  all  mesh  points.  For  components  not  evaluated  at  that 

mesh  point,  everything  on  the  right-hand  side  of  eq.  (10)  except  for"  S^e^ 

is   zero  and  in  that  term  S^  is  the  identity.   If  we  denote  the  vector  [e^, 

ez]   by  e  ,  we  get  a  difference  equation  for  e„  of  the  form 
n    J      n*    °  n 


where 


n+i     n    nnn  n 


-w 


■S 

0 

0 

Sn 

syy 

n 

fz 

gy 


,zz 


(11) 
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and  d^  is  a  sum  of  truncation  and  interpolation  errors.  The  latter  are 
bounded  by  °(h£  )  when  p-th  order  integration  and  (p-l)-st  order 
extrapolation  methods  are  used.  All  that  is  needed  is  to  prove  the 
stability  of  the  matrices  appearing  in  eq.  (11),  and  this  follows 
immediately  from  an  application  of  the  theory  in  {2}  since  the  matrices  S^ 
and  S^,  which  are  either  the  matrices  of  an  integration  method  or  identity 
matrices,  satisfy  the  stability  property  individually  and  the  matrix  S™  is 
bounded  because  its  four  components  are  bounded.  Stability  and  the  fact 
that  d^J  is  0(hP+1)  is  sufficient  to  prove  that  en  is  0(hp)  where  h  is  the 
maximum  stepsize. 

5.   VARIABLE-STEP  METHODS 

Very  few  problems  are  such  that  the  user  can  predetermine  the  stepsize 
to  be  used  throughout  an  integration,  so  most  modern  codes  vary  the  order 
and  stepsize  during  the  integration  based  on  local  error  estimates.  In 
this  section  we  look  at  the  problem  of  varying  the  stepsize  in  a  multirate 
method. 

There  are  problems  for  which  it  is  known  that  some  components  are 
always  faster  than  others,  and  for  these  problems  it  may  be  possible  to 
permanently  set  the  ratio  of  the  stepsizes.  Then  it  is  a  question  of 
choosing  one  step  size  and  letting  the  other  adjust  correspondingly.  A 
typical  automatic  code  performs  a  single  step  integration,  estimates  the 
error,  and  then  decides  whether  to  accept  the  step  because  the  estimated 
error  is  within  tolerance,  or  to  reject  it  and  repeat  it  with  a  smaller 
stepsize.  A  multirate  method  with  a  fixed  ratio  between  stepsizes  in 
different  components  could  be  viewed  as  a  type  of  cyclic  method  with  a 
cycle  equal  to  the  largest  stepsize  used  in  any  component.  Error 
estimations  could  be  made  over  that  cycle  and  either  the  cycle  accepted  or 
rejected  and  repeated.  Unfortunately  this  introduces  two  inefficiencies: 
the  loss  of  a  relatively  large  amount  of  work  when  a  cycle  is  rejected,  and 
the  need  to  back-up  a  number  of  steps  in  some  components.  This  means  that 
a  considerable  number  of  additional  past  values  must  be  saved. 

If  large-scale  back-up  and  loss  of  work  is  to  be  avoided,  errors  must 
be  estimated  as  each  step  is  made.   In  that  case,  there  does  not  seem  to  be 
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any  great  reason  for  fixing  the  stepsize  ratio;  indeed  in  general  we  cannot 
decide  a  priori  on  a  ratio  of  stepsizes,  some  components  may  be  the  fastest 
at  some  times  and  the  slowest  at  others.  Three  approaches  to  the 
organization  of  automatic  multirate  methods  have  been  tried.  In  all 
approaches  the  local  error  is  estimated  for  each  component  separately,  and 
then  the  step  in  that  component  is  either  accepted  or  rejected.  If  it  is 
rejected,  it  will  be  tried  with  a  smaller  stepsize  later.  The  approaches 
are: 

An  "event-driven"  model  in  which  each  mesh  point  is  viewed  as  an 
event,  and  the  next  event  in  time  is  chosen  for  execution. 

A  "recursive"  view  in  which  all  equations  are  attempted  with  the 
largest  possible  stepsize,  and  those  that  are  rejected  are 
integrated  with  a  pair  of  recursive  calls  at  half  the  stepsize. 

A  "largest-step-first"  approach  which  retains  some  of  the  advantages 
of  the  recursive  approach  but  reduces  function  evaluations. 

The  first  approach  was  tried  in  an  experimental  code  described  in  {5>. 
Each  component  of  the  system  was  integrated  with  a  possibly  different 
stepsize  and  order.  They  were  selected  by  the  integrator  for  each 
component  independently  of  the  behavior  of  the  other  components. 
Consequently,  at  any  given  time  each  component  i  had  been  integrated  to  a 
different  t  value  tj,  and  the  integrator  had  suggested  a  different  stepsize 
h^  for  the  next  step  for  each.  The  basic  strategy  was  to  select  the 
component  with  the  smallest  value  of  t^  +  hj  as  the  component  to  be 
integrated  next.  It  was  integrated  one  step,  its  error  estimated,  and  a 
value  for  its  next  stepsize  recommended.  The  idea  behind  this  is  that  the 
extrapolation  performed  in  the  other  components  is  to  points  within  the 
range  of  their  next  recommended  steps  and  within  this  range  the 
extrapolation  should  be  reasonably  accurate.  The  snag  to  this  arrangement 
is  that  the  recommended  step  may  be  too  large.  If  the  recommended  step  is 
far  too  large,  the  extrapolation  can  be  badly  in  error  causing  large  errors 
to  be  propogated  to  other  components.  The  difficulties  do  not  stop  there. 
If  a  recommended  step  turns  out  to  be  too  large  and  must  be  greatly 
reduced,  we  will  find  that  we  need  to  approximate  a  fast  component  at  a  t 
value  that  is  too  far  back  to  safely  extrapolate  backwards.   Consequently, 
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an  arbitrary  amount  of  back-up  may  be  necessary.  Clearly,  It  Is  not 
feasible  to  save  information  for  backing  up  to  arbitrary  earlier  values,  so 
we  tried  various  techniques  to  prevent  step  failures.  Consider  the 
integration  of  z  for  a  moment.  The  recommended  step  may  be  too  large 
because  of  a  change  of  behavior  in  g  that  is  not  predicted  by  earlier 
values  of  z,  or  by  a  change  in  the  behavior  of  y  that  couples  into  z 
through  g.  The  first  can  be  predicted  by  doing  a  trial  integration  of  z 
with  a  simple  extrapolation  for  y.  We  used  this  to  see  if  the  recommended 
step  seemed  reasonable.  It  is  clearly  expensive  because  it  almost  doubles 
the  number  of  integration  steps.  The  second  source  of  problems  was 
examined  by  trying  to  estimate  the  effect  of  changes  in  the  behavior  of  y 
on  z.  We  assumed  that  the  stepsize  chosen  for  z  had  been  based  on  the 
current  knowledge  of  y.  This  is  equivalent  to  saying  that  it  is  based  on 
the  extrapolated  value  of  y.  When  y  was  integrated,  the  predictor- 
corrector  difference  is  the  difference  between  the  extrapolated  value  of  y 
and  the  new  value.  The  effect  of  this  difference  on  the  z  integration  was 
estimated.  This  required  knowledge  of  the  Jacobian  g^.  The  method 
consisted  of  integrating  one  component  one  step  using  the  selection 
algorithm  given  above,  and  then  determining  if  the  recommended  steps  of  any 
other  components  had  to  be  reduced  because  of  the  error  in  the  step  just 
performed.  The  amount  of  work  was  made  manageable  for  a  system  of  several 
equations  by  keeping  the  sparsity  structure  of  the  Jacobian.  Only 
components  which  had  nonzero  entries  in  the  row  of  the  Jacobian 
corresponding  to  the  equation  just  integrated  had  to  be  considered.  The 
results  of  numerical  tests  of  this  technique  indicated  that  it  is  feasible 
if  the  evaluation  costs  of  g  are  sufficiently  high  and  the  ratio  of 
stepsizes  that  could  be  used  in  different  components  was  high  enough.  For 
a  number  of  test  problems  of  6  to  20  equations,  the  number  of  individual 
evaluations  of  derivatives  for  the  multirate  method  was  about  half  that  of 
the  standard  method.   However,  the  overhead  was  very  much  higher. 

The  principal  difficulty  in  the  event-driven  method  is  due  to  the 
integration  of  fast  components  before  one  is  certain  that  the  selected 
stepsize  for  the  slow  components  is  not  too  large.  Recall  from  Section  2 
that  we  said  the  coupling  between  fast  and  slow  components  is  necessarily 
small  and  it  might  make  sense  to  integrate  the  slower  component  first.   In 
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the  recursive  method  we  do  just  that.  All  components  are  Integrated  using 
the  largest  recommended  stepsize.  Then  the  error  in  each  is  estimated  and 
the  results  for  those  with  small  errors  are  accepted.  The  remainder  are 
reintegrated  using  two  steps  of  the  half  size  by  the  same  technique.  The 
advantages  of  the  recursive  approach  are  that  it  is  conceptually  simple  and 
not  too  difficult  to  program  even  if  recursion  is  not  provided  in  the 
language.  Also,  more  accurate  interpolation  rather  than  extrapolation  is 
used  to  get  the  intermediate  values  of  the  slow  variables  when  integrating 
the  faster  ones.  Its  disadvantage  is  that  there  are  many  unnecessary 
function  evaluations,  integrations,  and  interpolations  for  those  components 
whose  stepsizes  have  to  be  reduced  several  times.  In  fact,  there  are 
almost  twice  as  many  evaluations,  integrations,  and  interpolations  because 
if  a  step  of  H/2  is  used  when  the  largest  step  is  H,  the  number  of 
unsuccessful  steps  is  2  -  1  to  the  number  of  successful  ones  2  . 

The  third  technique  uses  the  same  principle  as  the  recursive 
technique,  that  is,  it  integrates  the  slowest  components  with  the  largest 
stepsizes  first.  A  maximum  stepsize  H  is  chosen,  and  all  stepsizes  are 
H/2  for  some  N.  Initially,  all  components  are  known  at  a  given  time  value 
tQ.  The  components  whose  stepsizes  are  H  are  integrated  first.  A  step  is 
halved  if  its  estimated  error  is  too  large.  Next  the  steps  of  size  H/2  are 
integrated,  and  so  on.  Finally,  the  step  with  the  smallest  size,  H/2  ,  are 
integrated.  This  process  is  repeated  for  all  components  whose  current  time 
value  is  least.  The  immediate  effect  of  this  is  to  integrate  the 
components  with  the  smallest  step  to  tQ  +  H/2'  .  Repeating  the  process 
again  causes  the  components  with  stepsizes  H/2^  "  '  and  H/2  to  be 
integrated,   in  that  order.   If  the  local  error  estimate  in  a  component  is 

M 

very  small,  the  stepsize  can  be  doubled.  However,  a  step  of  size  H/2  for 
M  >  0  can  only  be  doubled  when  the  time  status  of  its  differential  equation 
is  on  a  mesh  point  corresponding  to  steps  of  size  H/2^    '•   The  effect  of 

M 

this  is  that  components  using  stepsize  H/2  are  evaluated  only  on  the  mesh 
points  Tq  +  mH/2M  for  m  -  0,  1,  ...,  2M.  When  time  TQ  +  H  is  reached,  a 
new  maximum  stepsize  H  is  chosen  on  the  basis  of  the  error  estimates  for 
the  slowest  components. 

The  last  method  appears  to  be  the  most  efficient  on  the  basis   of 
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preliminary  tests.  Note  that  it  has  one  hidden  additional  cost.  It  is 
necessary  to  extrapolate  the  values  of  the  fast  components  in  order  to 
evaluate  the  derivatives  of  the  slower  components  which  are  integrated 
first.  Because  of  the  assumed  small  coupling  between  the  fast  and  slow 
components,  we  have  been  using  very  low  order  extrapolation  such  as  linear 
methods,  even  when  higher-order  integration  schemes  are  used. 

6.  STABILITY,  STIFFNESS,  AND  IMPLICIT  METHODS 

It  must  be  stressed  that  the  existence  of  slow  and  fast  components  has 
nothing  to  do  with  stiffness  per  se.  A  system  may  be  neutrally  stable  and 
still  profit  from  multirate  methods,  but  the  existence  of  stiffness  can  be 
a  complicating  factor  because  of  the  need  to  use  implicit  integration 
methods  and  to  solve  nonlinear  equations  at  each  step.  When  a  stiff 
problem  is  to  be  solved,  the  Jacobian  of  all  equations  being  integrated  to 
the  same  mesh  point  must  be  formed,  and  this  Jacobian  must  be  used  in  a 
quasi-Newton  iteration  for  the  solution  of  the  Implicit  integration 
formulas.  In  the  event  approach  above,  the  values  of  the  other  variables 
will  be  extrapolated  and  fixed.  In  the  recursive  method,  all  components 
with  larger  stepsizes  can  be  interpolated,  and  the  system  is  solved  for  all 
components  with  the  current  and  smaller  stepsizes.  In  the  largest-step- 
first  method,  values  of  variables  with  smaller  stepsizes  are  not  yet  known. 
Consequently  they  must  either  be  extrapolated — a  risky  business  in  stiff 
equations — or  can  be  held  constant.  This  corresponds  to  a  zero-th  order 
extrapolation  and  seems  to  be  the  best  choice.  In  systems  that  arise  from 
PDEs,  there  is  also  the  possibility  that  the  variable  with  the  smallest 
time  step  can  be  ignored  by  reverting  to  a  spatial  discretization  based  on 
a  coarser  mesh  or  finite  element  discretization.  This  is  related  to  the 
multigrid  techniques  discussed  in  {6}. 

7.  CONCLUSION 

The  use  of  multirate  methods  is  attractive  when  function  evaluation 
costs  are  very  high  or  if  there  is  a  high  degree  of  sparsity.  However,  the 
organization  of  an  automatic  code  is  not  simple,  and  it  is  very  easy  to 
allow  the  overhead  to  become  a  major  part  of  the  cost.  A  particular 
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problem  is  that  of  avoiding  back-up  over  several  steps  at  a  step  failure, 
and  the  counter-intuitive  approach  of  integrating  the  slowest  components 
first  seems  to  be  most  effective  in  this. 
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